%pylab inline
from pylab import *
import numpy as np # convert list to array
### Importing files from different folder
### https://stackoverflow.com/questions/4383571/importing-files-from-different-folder
import sys
sys.path.append('/home/cserti/programok/python/DOMAIN_COLORING')
from complex_func_plot_CsJ import * #notebook written by J. Cs. to plot complex functions
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import * #3D-s ábrák alcsomagja
from ipywidgets import * #interaktivitáshoz szükséges függvények
from matplotlib import cm
#from scipy import special
from scipy.special import jn,fresnel
from colorsys import hls_to_rgb
from matplotlib.colors import hsv_to_rgb
#import mpmath
The diffraction of waves by a half-plane, Sommerfeld's exact solution: ¶
The diffraction of a plane wave by a semiinfinite plane sheet is treated rigorously by Sommerfeld.
Original work: Sommerfeld, A 1896 ‘Theorie der diffraction’ Math. Ann. 47 317-374
$E_z$ polarizitation:¶
Here we assume that for the incident light $E_x=E_y=H_z = 0$, ie, the electric field is perpendicular to the incident plane called E-polarization. From now we introduce the notation $\psi(\mathbf{r}) = E_z(\mathbf{r})$.
The incident plane waves diffracting on an infinitely thin, half-plane defined by $x>0, y =0$ in polar coordinates ($x = \cos \varphi, y = \sin \varphi$, and $\varphi = [0,2\pi]$) is given by
$$ E_z(\mathbf{r}) \equiv \psi(\mathbf{r}) = e^{-ikr \cos\left(\varphi - \varphi_0\right)}, $$where the wave coming from the direction $\varphi_0$ to the origin and $k=2\pi/\lambda$ is the wave number.
Sommerfeld's famous result (may be found in Born & Wolf, 7th Ed., chapter XI, p. 633. see here),
the exact solution of this diffraction problem is
,
where $u_-(\mathbf{r}) = -\sqrt{2kr} \cos\frac{\varphi-\varphi_0}{2} $, $u_+(\mathbf{r}) = -\sqrt{2kr} \cos\frac{\varphi+\varphi_0}{2} $, and \begin{align} G(u) &= e^{-i u^2}\, \int_u^\infty e^{i t^2}\, dt = - e^{-i u^2}\sqrt{\frac{\pi }{2}} \left[ C\left(u \sqrt{\frac{2}{\pi}}\right) +i S\left(u \sqrt{\frac{2}{\pi}}\right)-\left(\frac{1}{2} +\frac{i}{2}\right) \right], \end{align} with $C,S$ are the Fresnel integrals defined according to python or Mathematica.
Boundary conditions:
$\alpha =1$ for Neumann, $\alpha =-1$ for Dirichlet boundary conditions, $\alpha =0$ for black screen.
The current distribution:¶
$$ \mathbf{j} = \mathrm{Im}\{\psi^* \mathrm{grad} \,\psi \} = \rho^2 \, \mathrm{grad}\, \chi , $$where $\psi (\mathbf r) = \rho(\mathbf r) e^{i \chi (\mathbf r)}$, with modulus $\rho$ and phase $\chi$. The current vector perpendicular to the contour lines of the phase $\chi (\mathbf{r})$.
The gradient of $\psi$ can be calculated analytically:
\begin{align} \frac{\partial \psi}{\partial x} &= - e^{-i\frac{\pi}{4}} \, \frac{e^{i k r}}{2r} \left\{ u_+ +\alpha \, u_- + i \, 2 k r \cos\varphi_0 \left[G(u_-) +\alpha \, G(u_+)\right] \right\}, \\ \frac{\partial \psi}{\partial y} &= -e^{-i\frac{\pi}{4}} \, \frac{e^{i k r}}{2r} \left\{ v_+ + \alpha \, v_- + i \, 2 k r \sin \varphi_0 \left[G(u_-) -\alpha \, G(u_+)\right] \right\}, \end{align}where $v_-(\mathbf{r}) = -\sqrt{2kr} \sin\frac{\varphi-\varphi_0}{2} $, $v_+(\mathbf{r}) = -\sqrt{2kr} \sin\frac{\varphi+\varphi_0}{2} $.
Note that from the Maxwell equations it follows that the magnetic field is perpendicular to $\mathrm{grad} \,\psi $, namely \begin{align} H_x & = \frac{1}{ik}\, \frac{\partial \psi}{\partial y}, \sim j_y \\ H_x & = -\frac{1}{ik}\, \frac{\partial \psi}{\partial x} \sim -j_x, \end{align} thus the magnetic field $\mathbf{H}$ is perpendicular to the current vector $\mathbf{j}$. Note that there could be a misprint in Eq. (28) on page 648 of Born & Wolf's book (7th Ed.), namely $H_y$ in the book should be multiply by a factor -1.
$H_z$ polarizitation:¶
The incident light $H_x=H_y=E_z = 0$, ie, the magnetic field is perpendicular to the incident plane called H-polarization. From now we introduce the notation $\psi(\mathbf{r}) = H_z(\mathbf{r})$.
Some other works:
In this paper the screen is rotated by $-90^\circ$ which we may denote by the rotation matrix $R$. Then the wave function in Berry's paper can be obtained from the Born-Wolf result $\psi (\mathbf{r}) $ shown above if we set $\varphi_0 = 3\pi/2$ and then rotate the wave function by $R$ according to the general rule: $\psi^{\left(\mathrm{Berry}\right)}\equiv \psi^\prime (\mathbf{r}) = \psi (R^{-1} \mathbf{r})$, i.e., $\varphi \to \varphi + 90^\circ$ and $r \to r$.
## The diffraction of waves by a half-plane, Sommerfeld's exact solution
def vszog(x,y):
# fi = [0,2*pi] normal definicio
return where(y<0,arctan2(y,x)+2*pi,arctan2(y,x))
## G(a) is defined in Eq. (25) on page 682 in Born-Wolf's book (7th editions)
def G_func(a):
tmp=fresnel(sqrt(2/pi)*a) # ---> FresnelS, FresnelC
res= -exp(-1j*a**2)*sqrt(pi/2)*(1j*tmp[0] + tmp[1] - (1+1j)/2)
return(res)
'''
Wave function:
$\psi$ defined above:
boundary ---> \alpha in the expression given above
boundary=0 # black screen boundary conditions
boundary=-1 # Dirichlet boundary conditions
boundary=1 # Neumann boundary conditions
'''
def psi_func_Born_Wolf(x, y, *params):
#E_z = psi is defined in Eq. (26) on page 682 in Born-Wolf's book (7th editions)
k=1 # wave number, r is in units of \lambda/(2 pi)
fi0, boundary = params # parameters for the function
r= sqrt(x**2+y**2) # sqrt(x**2+y**2)
fi=vszog(x,y)
tmp=sqrt(2*k*r)
um = -tmp*cos((fi-fi0)/2)
up = -tmp*cos((fi+fi0)/2)
psi=exp(-1j*pi/4)/sqrt(pi)*exp(1j*k*r)*(G_func(um)+boundary*G_func(up))
return(psi)
def psi_func_Born_Wolf_Hz(x, y, *params):
#H_z = psi is defined in Eq. (45) on page 687 in Born-Wolf's book (7th editions)
k=1 # wave number, r is in units of \lambda/(2 pi)
fi0, boundary = params # parameters for the function
r= sqrt(x**2+y**2) # sqrt(x**2+y**2)
fi=vszog(x,y)
tmp=sqrt(2*k*r)
um = -tmp*cos((fi-fi0)/2)
up = -tmp*cos((fi+fi0)/2)
psi=exp(1j*pi/4)/sqrt(pi)*exp(1j*k*r)*(G_func(um)-boundary*G_func(up))
return(psi)
def psi_func_Born_Wolf_current(func, x,y,*params):
fi0, boundary = params # parameters for the function
k=1 # wave number, r is in units of \lambda/(2 pi)
r=sqrt(x**2+y**2)
fi=vszog(x,y)
#fi=arccos(x/r)
#fi=arctan2(y,x)
psi = func(x,y,*params)
tmp=sqrt(2*k*r)
um = -tmp*cos((fi-fi0)/2)
up = -tmp*cos((fi+fi0)/2)
vm = -tmp*sin((fi-fi0)/2)
vp = -tmp*sin((fi+fi0)/2)
fact = -exp(-1j*pi/4)*exp(1j*k*r)/(2*r)
tmp1 = up + boundary*um
tmp2 = 1j* 2*k * r* cos(fi0)* (G_func(um)+boundary*G_func(up))
gradpsi_x = fact*(tmp1+tmp2)
tmp1 = vp + boundary*vm
tmp2 = 1j* 2*k * r* sin(fi0)* (G_func(um)-boundary*G_func(up))
gradpsi_y = fact*(tmp1+tmp2)
jx= imag(conjugate(psi)*gradpsi_x)
jy= imag(conjugate(psi)*gradpsi_y)
# for test
#jx=-y
#jy=x
return(jx,jy)
'''
arctan2(y,x) does not give the right angle, according to the above definition of \varphi:
'''
rp=[[1,1],[-1,1],[-1,-1],[1,-1]] # points in the four quadrant
szog=[]
[szog.append(180/pi*arctan2(rp[i][1],rp[i][0])) for i in range(len(rp))]
rp,szog
# Now it is ok
rp, [180/pi*vszog(rp[i][0],rp[i][1]) for i in range(len(rp))], 180/pi*vszog(0,-1)
G_func(1+2j)
x,y =(-3,2)
for i in array([-1,0,1]):
print(i,psi_func_Born_Wolf(x, y, *[pi/2,-1]))
x,y =(-3,2)
for i in array([-1,0,1]):
print(i,psi_func_Born_Wolf_current(psi_func_Born_Wolf,x, y, *[pi/2,-1]))
x,y =(-3,2)
for i in array([-1,0,1]):
print(i,psi_func_Born_Wolf_Hz(x, y, *[pi/2,-1]))
x,y =(-3,2)
for i in array([-1,0,1]):
print(i,psi_func_Born_Wolf_current(psi_func_Born_Wolf_Hz,x, y, *[pi/2,-1]))
$ \left\{-i k \cos (\text{fi0}) G(\text{uum}) e^{i k r}-i k \sigma \cos (\text{fi0}) G(\text{uup}) e^{i k r}+\frac{\sigma e^{i k r} \sqrt{k r} \cos \left(\frac{\text{fi0}-\var phi }{2}\right)}{\sqrt{2} r}+\frac{e^{i k r} \sqrt{k r} \cos \left(\frac{\text{fi0}+\var phi }{2}\right)}{\sqrt{2} r},-i k \sin (\text{fi0}) G(\text{uum}) e^{i k r}+i k \sigma \sin (\text{fi0}) G(\text{uup}) e^{i k r}-\frac{\sigma e^{i k r} \sqrt{k r} \sin \left(\frac{\text{fi0}-\var phi }{2}\right)}{\sqrt{2} r}+\frac{e^{i k r} \sqrt{k r} \sin \left(\frac{\text{fi0}+\var phi }{2}\right)}{\sqrt{2} r}\right\} $
def test_fn(x,y,*params):
z=x+1j*y
res = z # for test
return(res)
Here we use the conventional color coding¶
See the color wheel based on a colormap using Python/Matplotlib
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
domain = [-3, 3, -3, 3]
Np=100
params = ''
title='$f(z)=z$'
#title=''
#fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
#fig, ax =plt.subplots(1,1, figsize=(5, 4), sharey=True)
fig, ax =plt.subplots(1,1, figsize=(5, 5), sharey=True)
wave_func_phase_color(test_fn, domain, Np, title, *params)
#wave_func_2D(test_fn, domain, Np, title, *params)
#ax.set_title(title);
wave_func_phase_contour(test_fn, domain, Np, title, *params);
#plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
show()
$E_z$ polarization¶
Dirichlet boundary conditions, boundary=-1¶
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [1.*pi/2, -1] # fi0, boundary
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
size = 2*pi
domain=array([-size,size,-size,size])
Np = 200
N_level = 15
title = '$E_z(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#title=''
#fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
#fig, ax =plt.subplots(1,1, figsize=(6, 5))
fig, ax =plt.subplots(1,1, figsize=(10, 10), sharey=True)
wave_func_2D(psi_func_Born_Wolf, domain, Np, title, *params)
#ax.set_title(cim)
wave_func_2D_contour(psi_func_Born_Wolf, domain, Np, N_level, title, *params);
#
plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [pi/2, -1] # fi0, boundary
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
size = 5*pi
domain=array([-size,size,-size,size])
Np = 200
title = '$E_z(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#title=''
#fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
fig = plt.figure(figsize=(12,6))
plt.subplots_adjust(wspace = 0.33 )
ax1=plt.subplot(1,2,1)
wave_func_2D(psi_func_Born_Wolf, domain, Np, title, *params)
#ax.set_title(cim)
plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
#ax1.axhline(y=0, xmin=0, xmax=domain[1], color='k', lw=7, ls='solid')
ym= -1*pi
ax1.axhline(y=ym, xmin=domain[0], xmax=domain[1],color='r', ls='--')
ax2=plt.subplot(1,2,2)
x0=6*pi
xv=linspace(-x0,x0,Np)
plot(xv,abs(psi_func_Born_Wolf(xv, ym, *params)))
ax2.axhline(y=1, xmin=xv[0], xmax=xv[-1],color='k', ls='--')
ax2.axhline(y=0.5, xmin=xv[0], xmax=xv[-1],color='k', ls='--')
ax2.set_title('$\Psi(\mathbf{r}), $'+ ' boundary = ' + str(params[1]))
ax2.set_xlabel(r'$x\, (\lambda/2 \pi)$', fontsize=15)
ax2.set_ylabel(r'${|\psi(x)|}^2$', fontsize=15)
#ax2.set_aspect('equal');
ax2.grid()
fig.tight_layout();
print('ym (in units of wave length)= ',ym/(2*pi))
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [1.*pi/2, -1] # fi0, boundary
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
size = 2.0*pi
domain=array([-size,size,-size,size])
Np = 200
title = '$\chi(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#title=''
fig, ax =plt.subplots(1,1, figsize=(8, 8), sharey=True)
wave_func_phase_color(psi_func_Born_Wolf, domain, Np, title, *params);
wave_func_phase_contour(psi_func_Born_Wolf, domain, Np, title, *params);
plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [1.*pi/2, -1] # fi0, boundary
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
size = 2.0*pi
domain=array([-size,size,-size,size])
Np=100
title = '$\chi(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#title=''
fig, ax =plt.subplots(1,1, figsize=(8, 8), sharey=True)
wave_func_current(psi_func_Born_Wolf, psi_func_Born_Wolf_current, domain, Np, title, *params);
#title(title)
wave_func_phase_contour(psi_func_Born_Wolf, domain, Np, title, *params);
plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [1.*pi/2, -1] # fi0, boundary
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
size = 2.0*pi
domain=array([-size,size,-size,size])
Np=100
title = '$\mathbf{j}(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#title=''
fig, ax =plt.subplots(1,1, figsize=(8, 8), sharey=True)
wave_func_current(psi_func_Born_Wolf, psi_func_Born_Wolf_current, domain, Np, title, *params);
#title(title)
plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [0.5*pi/2, -1] # fi0, boundary
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
size = 3*pi
domain=array([-size,size,-size,size])
Np=100
N_level = 15
title = '$E_z(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#title=''
#fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
#fig, ax =plt.subplots(1,1, figsize=(6, 5))
fig, ax =plt.subplots(1,1, figsize=(10, 9), sharey=True)
#plot_domain(domaincol_co, psi_func_Born_Wolf, domain, Np, title, *params, daxis=True);
wave_func_phase_color(psi_func_Born_Wolf, domain, Np, title, *params)
wave_func_phase_contour(psi_func_Born_Wolf, domain, Np, title, *params);
#ax.set_title(cim)
#wave_func_current(psi_func_Born_Wolf, psi_func_Born_Wolf_current, domain, Np, title, *params);
#
plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [1.*pi/2, -1] # fi0, boundary
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
size = 5*pi
domain=array([-size,size,-size,size])
Np=100
title = '$E_z(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
#fig, ax =plt.subplots(1,1, figsize=(5, 4), sharey=True)
#figsize=(12,8)
wave_func_3D(psi_func_Born_Wolf, domain, Np, title, *params)
#ax.title(title)
#plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
#plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
show()
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [1.0*pi/2, -1] # fi0, boundary
size = 5*pi
domain=array([-size,size,-size,size])
Np=100
title = '$\chi(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#title=''
figure(figsize=(12,8))
plt.subplot(2, 2, 1)
plot_domain(domaincol_c, psi_func_Born_Wolf, domain, Np, title, *params, daxis=True)
plt.subplot(2, 2, 2)
plot_domain(domaincol_m, psi_func_Born_Wolf, domain, Np, title, *params, daxis=True)
plt.subplot(2, 2, 3)
plot_domain(domaincol_co, psi_func_Born_Wolf, domain, Np, title, *params, daxis=True)
plt.tight_layout();
Neumann boundary conditions, boundary=1¶
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [pi/2, 1] # fi0, boundary
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
size = 3*pi
domain=array([-size,size,-size,size])
Np = 200
N_level = 15
title = '$E_z(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#title=''
#fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
#fig, ax =plt.subplots(1,1, figsize=(6, 5))
fig, ax =plt.subplots(1,1, figsize=(10, 10), sharey=True)
wave_func_2D(psi_func_Born_Wolf, domain, Np, title, *params)
#ax.set_title(cim)
wave_func_2D_contour(psi_func_Born_Wolf, domain, Np, N_level, title, *params);
#
plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
Black screen boundary conditions, boundary=0¶
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [pi/2, 0] # fi0, boundary
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
size = 2*pi
domain=array([-size,size,-size,size])
Np = 200
N_level = 15
title = '$E_z(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#title=''
#fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
#fig, ax =plt.subplots(1,1, figsize=(6, 5))
fig, ax =plt.subplots(1,1, figsize=(10, 10), sharey=True)
wave_func_2D(psi_func_Born_Wolf, domain, Np, title, *params)
#ax.set_title(cim)
wave_func_2D_contour(psi_func_Born_Wolf, domain, Np, N_level, title, *params);
#
plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
$H_z$ polarization¶
### half-plane scattering, Sommerfeld solution, Born-Wolf
params = [pi/2, -1] # fi0, boundary
# domain --> domain to be plotted:
# x --> (domain[0],domain[1]),
# y --> (domain[2],domain[3])
size = 2*pi
domain=array([-size,size,-size,size])
Np = 200
N_level = 15
title = '$H_z(\mathbf{r}), $'+ ' boundary = ' + str(params[1]) + r'$, \,\, \varphi_0 = $' + str(round(params[0]*180/pi,1))
#title=''
#fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
#fig, ax =plt.subplots(1,1, figsize=(6, 5))
fig, ax =plt.subplots(1,1, figsize=(10, 10), sharey=True)
wave_func_2D(psi_func_Born_Wolf_Hz, domain, Np, title, *params)
#ax.set_title(cim)
wave_func_2D_contour(psi_func_Born_Wolf_Hz, domain, Np, N_level, title, *params);
#
plt.hlines(y=0,xmin=0, xmax=domain[1], colors='k', lw=7, linestyles='solid');
